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Abstract. We define the general Hill system and briefly analyze its dynamical behavior. A 
particular Hill system representing the interaction of a Keplerian binary system with a normally 
incident circularly polarized gravitational wave is discussed in detail. In this case, we compute 
the Poincare-Melnikov function explicitly and determine its zeros. Moreover, we provide numerical 
evidence in favor of chaos in this system. The partially averaged equations for the Hill system are 
used to predict the regular behavior of the Keplerian orbit at resonance with the external radiation. 
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1 Introduction 



In a previous paper on gravitational ionization [[!]], we discussed — among other things - 
the long-term dynamical evolution of a Keplerian binary system perturbed by a normally 
incident circularly polarized gravitational wave. If the external radiation is approximated by 
a monochromatic plane wave, then the dynamical system representing the relative motion 
of the binary is very similar to the well-known Hill system of celestial mechanics. Moreover, 
in our numerical work on such a system, there appeared some preliminary evidence in favor 
of Hamiltonian chaos (see figure 1 of The purpose of the present work is to extend the 
notion of a Hill system to include the external perturbation caused by a gravitational wave. 
The generalized Hill system — which we refer to simply as "the Hill system" due to its close 
resemblance to the one originally introduced by Hill in his researches on the lunar theory 
- is presented in section |2|. For the sake of simplicity, we then choose a particular case for 
detailed analysis and discuss its dynamics near resonance in section |3|. Moreover, in section |] 
the relevant Poincare-Melnikov function is computed in this case and its countable infinity 
of zeros are studied. The possible relevance of this infinity of simple zeros to the existence 
of chaos in the Hill system is discussed. Numerical evidence for such chaos is presented in 
section [5] and a specific chaotic orbit is described in some detail. In section ^ we turn our 
attention to the physics of such orbits in the inertial frame in connection with the possible 
astrophysical relevance of our results. Section ^ is a discussion. Calculational details are 
relegated to the appendices. 



Consider a Keplerian binary system under the influence of an external gravitational per- 
turbation. In fact, no binary system in the universe is totally isolated as a consequence of 
the universality of the gravitational interaction. The binary is generally affected by other 
masses as well as gravitational radiation. We are interested in a particular form of this 
interaction that corresponds to Hill's celebrated contribution to the theory of the motion of 
the Moon @. 

Let us therefore assume that the predominant effect of the external perturbation is a 
linear tidal force in the equation of relative motion, so that this equation can be written as 



where r = |X|, k Q is a positive constant and i,j = 1,2,3. Here the tidal matrix K(t) is in 
general symmetric and traceless. Moreover, we assume that K 13 = K 2 z = 0; indeed, a tidal 
matrix of this form is consistent with our further requirement that the perturbed relative 
motion be confined to the (X, F)-plane. To arrive at the Hill system, we need to restrict the 
form of K(t) even further; to this end, let us suppose the following general form for K: 



2 The Hill System 



al 2 X l LI* 



K^X* 




alt 2 



K\\ = £o + £+ cos fit — sin fit, 
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K12 = C+ sin Qt + £_ cos fit, 
K22 = Co ~ C+ cos ^ + C- s i n ^? 

^33 = -2&, (2.2) 

where fi, Co and C± are real constants. 

It is possible to transform the equation of relative motion ( |2.1|) to a rotating coordinate 
system via X 1 = S^x- 7 , where 5* is an orthogonal matrix given by 



S 



cos (-fit) -sin (-fit) 
sin (ffit) cos(|fit) 
1 



(2.3) 



that represents a rotation by an angle |fit about the Z-axis. The equation of relative motion 
with respect to the rotating frame of reference with coordinates x = (x, y, z) is then given 
by the autonomous Hill system 

d 2 x dy 1 2 k x , , 

^- n Tt~l nx + — = "(^ + ^), 

§ + ^-\^ + ^ = _ {k2lX + k22yl (2 . 4) 

and z = 0, where r = |x| and A; = S^^^KS is a constant matrix given by k\\ = Co + £+> 
k 12 = ^22 = Co ~ C+j ^13 — ^23 = and k 33 = — 2Co- It is important to note that the sign 
of fi has not been specified; in fact, fi can be positive or negative. In Hill's original system, 
Co = — fi 2 /8, C+ = — 3fi 2 /8 and £_ = 0, corresponding to the circular motion of the Earth- 
Moon system about the Sun with frequency fi/2. Hill's variational orbit, which is a periodic 
solution of the Hill system, has played a major role in the lunar theory [|, |[ [5], ||. In addition 
to the original system envisaged by Hill, the general case includes a Kepler ian binary system 
that is tidally perturbed by a normally incident circularly polarized gravitational wave M. 

The Hill system may be obtained from the Hamiltonian H(p x ,p y , x, y), 

H = - p 2 -h- Iql + Ijk-zV', (2.5) 
2 r 2 2 3 

where the canonical momenta are given by p x = x — Qy/2, p y = y + Qx/2 and L z = xp y —yp x . 
Here x = dx/dt and L z is the component of relative orbital angular momentum normal to 
the orbital plane. It is interesting to introduce in this plane polar coordinates (r, 6) such 
that x = rcos6 and y = rsin^. Then with respect to the canonical variables (p r ,pe,r,6), 
equation ( [2.5|) can be written as 

H = \{£ + 1 %- 1 J- \toP» + \* (£0 + C + cos 29 + C- sin 26) , (2.6) 

where p r = x • p/r and p# = L 2 . In the absence of the external perturbation, i.e. when fi, 
Co and C± vanish, the relative motion of the binary can be described in terms of a Keplerian 
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ellipse. In the presence of the perturbation, the relative motion in the rotating system of 
reference can be described at each instant of time in terms of an osculating ellipse. That 
is, if the perturbation is turned off at an instant of time t, the subsequent motion is along 
the osculating ellipse at t. The energy E < and orbital angular momentum pg > of 
the osculating ellipse fix its semimajor axis a = —ko/(2E) and eccentricity e, < e < 1, 
Pg = a(l— e 2 ). The position on the osculating ellipse at time t is measured from the periastron 
by the true anomaly v , while the orientation of this whole osculating ellipse is determined in 
the orbital plane by an angle g such that 6 = v + g. The equation of the osculating ellipse 
is then given by r = pg/(l + ecosv) = a(l — ecosu), where u is the eccentric anomaly 
Moreover, p r pg = esinv, (1 — ecosw) sinv = (1 — e 2 )a sintt and the mean anomaly £ is given 
by £ = u — esinw. Only positive square roots are considered throughout this paper. The 
parameters of the osculating ellipse are closely related to Delaunay's action-angle variables. 
In fact, the Delaunay elements are defined by L = (koa)^, G = pg, £ and g for noncircular 
orbits. It is proved in || that the change of variables (p r ,Pe, r, 0) — > (L, G, £, g) is canonical. 
In terms of these planar Delaunay elements, we then have the Hamiltonian for the Hill 
system in the form 



H(L,G,£,g) 



fc ° - -QG 



2L 2 



(2.7) 



Here Q = r 2 is given by 



Q{L,G,£,g) = a 2 



, ,3 2 ~ COSI/* j . , 

! + o e - 4 E— —Juiye) 



(2.8) 



and C = r 2 cos 28 and S = r 2 sin 26 are given in Delaunay's elements by 



C + iS = a 2 exp(2ig) 



-e 2 + £ - fc exp(iu£) + K v _ exp(-iui 



(2.9) 



where 



and 



K v ± =~v{A v ±B v 



A v (e) 
BJe) 



v 2 e 2 



v 2 e 2 



[2ue(l - e 2 )J' v (ue) - (2 - e 2 )J v {ve% 
(\-e 2 )^{eJ' v (ve)-v(\-e 2 )J v (ye)\. 



(2.10) 
(2.11) 



It is interesting to note that A v = A_ u and B v 



-B_ u , so that Kl u 



It follows from our previous work |T|] that the Kolmogorov- Arnold- Moser ( "KAM" ) theory 
is applicable to the general Hill system and that for sufficiently small £o and £± the motion 
is always bounded. Moreover, our previous work (cf. figure 1 in 0) indicates the presence 
of Hamiltonian chaos under certain circumstances. The purpose of the present paper is to 
investigate further the nature of this chaotic behavior. 
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To simplify the analysis, we will assume in the following that £o = and that £ + = 
eafl 2 cos ipo and £_ = eaQ 2 siny? - For Q > 0, this choice of parameters corresponds to the 
tidal perturbation produced by an incident right circularly polarized (i.e. positive helicity) 
plane monochromatic gravitational wave of amplitude ate, frequency Q and phase constant 
ipo propagating along the z-axis. On the other hand, we can let Q — > — Q in equations 
(2)-(7) and the resulting system with the same choice of parameters would correspond to 
left circularly polarized (i.e. negative helicity) radiation of frequency Q > 0; in this way, our 
analysis covers both cases of circular polarization. The transverse nature of gravitational 
radiation implies that for this case of normal incidence the binary motion remains planar. 
The deviation of the curved spacetime metric in the presence of gravitational waves from the 
flat Minkowski metric is characterized by the perturbation parameter e, < e < 1. Efforts 
are under way in several laboratories to detect gravitational waves from astrophysical sources 
with an expected amplitude of e ~ 10~ 20 . Moreover, the circular polarization amplitude a 
is such that |a| ~ 1. To simplify matters even further, we shall set <^ = 0. The resulting 
Hill system has been discussed at length in our previous papers [7[ ^ . More generally, we 
have investigated the long-term nonlinear evolution of a Keplerian binary system perturbed 
by external long-wavelength gravitational waves as well as internal gravitational radiation 
damping []7|, |9], [10]. In fact, gravitational ionization provided the original motivation for 



our analysis [lTJ. The term "ionization" is derived from a Greek word meaning "to go"; 



therefore, the concept of ionization need not be restricted to electrically charged systems 
such as in atomic physics. We have shown that in the case under consideration ionization 
does not in fact occur for e < ek am - This is an important consequence of the KAM theory 
and implies, on the physical side, that the circularly polarized gravitational wave does not 
steadily deposit energy into the orbit. Indeed, in the interaction of gravitational waves with a 
binary system the energy in general flows back and forth between the waves and the binary. 
Further discussion of gravitational ionization is contained in section |B| in connection with 
certain large-scale chaotic behavior of the binary orbit for e > e^, where e c ^ corresponds to 
Chirikov's resonance-overlap condition. 



3 Isoenergetic Reduction and Dynamics Near Reso- 
nance 

With the simplifications already mentioned, the equations of motion derived from the Hamil- 
tonian (|2.7| ) are given by 



L = 


1 o2 dC 

~2 ean TV 




G = 


1 o2 dC 




t = 


k 2 1 o2 9C 




9 = 


L 1 o2 <9C 

-2 n+ r an dG- 


(3.1 
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These equations provide a classic example of a perturbation problem in mechanics. The 
unperturbed system is expressed in action-angle variables and is completely integrable. We 
are interested in the dynamics of the perturbed system for small e. This is the "fundamental 
problem in dynamical systems", according to Poincare H]. 



System ( |3.1| ) is a 2-degree-of-freedom autonomous Hamiltonian system. Thus by a well- 
known method that we will now describe, it can be reduced to a l|-degree-of-freedom system 
once the total energy of the system is fixed. Let us describe the reduction technique for the 
general case. In fact, we will consider the system 

dH 
OH 

p = —^-(q,p,$,i), 

<q,P,4,I), (3-2) 



where the Hamiltonian function H is assumed to be periodic in the angular variable and 
there is some region TZ of the phase space in which the function dH/dl does not vanish. 
Under these assumptions, we can solve for / as a function of the remaining variables on an 
energy surface, 

{{q,p,#,I):H(q,p,&,I) = h}, 

to obtain 

I = K(q 7 p^,h). 

Moreover, because dH/dl does not vanish in TZ, the variable is either increasing or de- 
creasing along all orbits in TZ. Thus, behaves like a temporal variable and can be taken to 
be a new independent variable. Indeed, the system 

dq _ dH dH dp __ dH dH 
is not singular in TZ. If •& i— > («(#), $(•&)) is a solution of the system (|3.3|), then 

PjTT 

= —(u0&), 5(0), 0, k(uO&),vO&), 0, h)). 
ol 

The solution 1 1— > 0(t) of this scalar differential equation can be used to obtain a solution of 
the original system (|3.2| ); namely, 



g(t)=w(0(t)), p{t)=v(&(t)), I(t) = k(q(t),p(t),#(t),h). 



To obtain a simpler form for the system ( |3.3| ), let us use the equation 

H(q,p,#,k(q,p,#,h)) = h 
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to obtain the identities 



dH OH OK _ dH dH dK _ 
dq dl dq dp dl dp 



and, in turn, system ( |3.3|) in the form 



d-d 9p W ' p ' ' h d$ dq U>P> v > n )- V A > 



Under our assumptions, the system (p.4|) is periodic in its independent variable. Moreover, 
system (3.4 ) is in the form of a "time-dependent" l^-degree-of- freedom Hamiltonian system 
with Hamiltonian —K(q,p, h). Let us now employ this reduction principle in the study of 
the system (|3.1|). 



For the Hamiltonian ( |2.7| ), we recall that if e is sufficiently small, then the KAM Theorem 
applies and all orbits in a region 1Z of an energy surface remain bounded. In particular, if TZ 
is a closed subset of a bounded region, then the function dC/dG is bounded over TZ. Thus, 
if e is sufficiently small, then g < along all orbits starting in TZ. Also, by the Implicit 
Function Theorem, if the constant energy surface is given by 

H(L,G,£,g) = h, 

then there is an implicit solution G = F(L,£,g) such that 

H(L,F(L,e,g),£,g) = h. 

Using the reduction principle given above, the reduced system in our case is given by 

M 8F (T 9 \ dL 9F (T 9 \ (**\ 

Tg = -dL {LJ > 9) > Tg=M^^ g) - (3 ' 5) 



To solve for G in the equation H(L, G, £, g) = h, let us suppose that in the corresponding 
equation ( p.7|) we have G = G + eG\ + 0(e 2 ) and then equate coefficients to obtain the 
solution up to first order in the small parameter. Let us recall here that the orientation 
of the background inertial coordinate system is so chosen that G > by convention. Our 
computations yield the values 

Go = ~ (h + ^ ) , Gi = ^ C(L, G , t,g)-~hi, 

where h and hi are constants such that h = h + ehx\ that is, h is the unperturbed energy 
in the rotating frame. Thus, we have that 

fg = - «*(§E<*. <*» *> * + °* <> + °^ 

^ = ean^(L, Go, £, g) + 0(e 2 ). (3.6) 
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Also, let us note that in view of equation ( |2.9| ) the l|-degree-of-freedom system ( |3.(j| ) is 
periodic with period ir relative to its independent variable. 



The unperturbed system ( |3.6| ) is a priori stable; that is, it has no hyperbolic orbits. In 
particular, there are no homoclinic or heteroclinic orbits that can be used to locate transverse 
homoclinic or heteroclinic points, and, in turn, chaotic invariant sets, after perturbation. 
Rather, if such transverse homoclinic or heteroclinic points exist in the perturbed system, 
they must be created directly from the perturbation. As is well known, this fact precludes a 
direct application of the usual "Melnikov" approach to determine the existence of transverse 
homoclinic or heteroclinic points. The problem is that the separatrix splitting distance 
is easily computed to be proportional to eM + 0(e 2 ), where M is the Melnikov integral. 
In the usual case, M is independent of e and thus the first term of the indicated series 
is the dominant term. However, in the a priori stable case, M depends on e and is in 
fact exponentially small. Therefore, it is not clear if eM(e) is the dominant term in the 
expression for the splitting distance. On the other hand, some rigorous results do exist for 
systems similar to those that are encountered in the Hill system. Consider, for instance, the 
rapidly forced pendulum given in the form 

(p + sin (p = e p sin(f2t/e), 

where e is a small parameter and p is a positive integer. If p > 3, then the Melnikov integral 
is exponentially small, but the first-order term is still dominant, see []T2 and [ |T3|j . 



While the rigorous analysis of the separatrix splitting problem remains unclear at this 
writing, we will carry through the first-order analysis for Hill's system. Once the rigorous 
analysis is settled at some future date, we hope that the formulation reported here will prove 
to be valuable. 

To proceed with the perturbation theory for Hill's system, we will use a standard tech- 
nique: we will "blow up" a resonance and partially average the resulting equations to produce 
a system, obtained by a change of coordinates, that has homoclinic orbits. We will then com- 
pute the Melnikov integral for this system. For a regular perturbation problem where the 
Melnikov function does not depend on the perturbation parameter, the simple zeros of this 
function would indicate the presence of chaos in the system. After a blow up at a resonance 
and a reparametrization of the system to slow time, the Melnikov function does depend on 
the perturbation parameter. It is for this reason that we are not able to make a rigorous 
argument that the existence of simple zeros implies that the system is chaotic. 



For system ( |3.6| ) that is periodic in g with period tt, resonance would occur when this 
period is commensurate with the period of "fast" motion in £, i.e. irQL 3 /k^. Therefore, the 
resonances are given by relations of the form 

m— = nil 

where m and n are relatively prime integers. This resonance condition fixes the magnitude 
of L; therefore, let L denote the resonant value of L at the (m : n) resonance. To blow up 



the resonance, we will use the change of coordinates given by 

2n 

L = L + e 1/2 p, £= 9 + V- 

m 

System ( |3.6| ), in the new coordinates, has the form 

d -f- = e^and + eatlp{C lL + -C lG ) + 0(e 3 / 2 ), 
dg m 

where the function C and its partial derivatives are evaluated at 

(Lo,G {L ),r] - 2ng/m,g). 

For notational convenience, let us define [i = e 1 / 2 and express the second-order approximation 
of system (|3.7|) in the form 

dg 

JL = p C p- p 2 (3(p,r],g), (3.8) 

where c is a constant given by c = 6n/ (mL ) and 

2n 2c 2n 

7{n,9)=<*®Ct, K(r),g) = aQ(C eL ^ C iG ), @(p, r/, g) = —p 2 + q>VL{Cl H C G ). 

m L m 

This system is in the correct form for temporal averaging. To this end, we define (7) to 
be the average of 7 over the periodic temporal variable g with period mir and A to be the 
deviation of 7 from its mean, 



(i)(v) ■ = — / i{v,g) d 9, 

mir Jo 

A(77,<7): = T (77,<7)-<7>(»7). (3.9) 

Moreover, let A denote the solution of the partial differential equation dA/dg = X with 
the side condition that its average should vanish; that is, (A) (77) = 0. Using the averaging 
transformation 

p = p + fiA(fj,g), rj = fj, 

system ( |3.8|) takes the form 

^ = v(i)(v) + pfpWv, g) - 4^) + o(/i 3 ), 

dg on 

p- = pcp + p 2 (cA(n, g) - p(j>, rj,g)) + 0(p 3 ). (3.10) 
dg 

It is useful to introduce dimensionless quantities J and F such that 

p = L J, A = L Q T] 
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then, the second-order approximation of system ( P-1C| ) has the convenient form 

dg Lq v 

JL = 6— fiJ + )i 2 (Q—T(r),g) - (3(L J,rj, g)). 
do "» m ' 



6n dT\ 
m dr)' 



m 



m 



(3.11) 



Recall formula ( |2.9|) and the fact that 

7 = aQCe(L ,G (L ), g + Vio), 

v m ' 

and note that (7) (77) is the real part of 

\ rmm , 2n 2n \ 

/ att(C £ {L ,G {L ), g + r),g)+iS e (L ,G (L ), g + r),g))dg 

mir Jo v m m ' 

I exp(2ig) [KK expUun — i — vg) — K v _ exp(— + i — vg) ) dg. 
Jo ~y v m m ' 



aVLa 2 



ran 



This integral vanishes unless n = 1 and m = u, in which case its value is iaQa 2 K™ exp(imr)). 
Using the fact that a = L 2 /ko, L = L at resonance and f2 = mk 2 LQ 3 , the average of 7 at 
the (m : 1) resonance can be written in the form 



Lo 



(l)(v) = —arnK™ sin mr). 



(3.12) 



If in system ( [3.11 ) we substitute using formula ( [3.12 ), change to the new angular variable 
a := mr) — tt and the slow temporal variable r := 6/zg, then we obtain the equivalent system 



da 



J + )i 



dr 

dJ am 



a + 7T T 



dr 6 



m ' 6/j,' 6 
K™ sin a + \i J 



^{LoJ, 



a + IT T 

m ' 6/i 



1 / (X + 7T T \ 8T rO + 7T T 

V m. Pill.' Fin V 
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m 6/i' da^ m 6/i 



(3.13) 
(3.14) 



Let us note that the unperturbed form of the system ( p,13| )- (|3.14| ), given by 

d<j 



dr 



t ^0 r • 

J , — = -dsma , 
dr 



is the phase-plane system for a mathematical pendulum. Here 

am 



6 



-K m 



(3.15) 



m = 1, 2, 3, . . ., is the order of the resonance, the constant a could be positive or negative 
and K™, 2K™ = m(A m + B rn ) as defined by equations ( |2.1(J| ) and ( |2.11| ), is a function of 
the eccentricity at resonance (eo). Some properties of K±(e) near e = and e = 1 are 
given in Appendix [A[ A graph of the functions K± for m = 1,2,3 is given in figure 1, 
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Figure 1: Plot of K™ for m = 1, 2, 3. Note that K™ and K™ meet at e = 1. Some of the 
properties of K± are discussed in Appendix [A]. 

which illustrates the analytical results of Appendix [A[ Our explicit calculations refer to the 
case of incident positive helicity radiation and hence only K™ is involved in the definition 
of 5 (see section |6|). It follows from inspection of figure 1 that depending upon the value of 
the eccentricity eo, 5 could be positive, negative or zero. If 5 > 0, then the phase portrait 
has hyperbolic saddle points at ( J , cr ) = (0, ±7r) that are connected by the heteroclinic 
solutions 

0" o (t) = 2 arcsin(tanh(5 1 ^ 2 r)) = 2 arctan(sinh(5 1 / 2 r)), 
J (r) = ±25 1/2 sech(5 1/2 r). 

If 5 < 0, then the phase portrait has hyperbolic saddle points at (Jo,cr ) = (0,0) and 
( Jq, Co) = (0, 27r) that are connected by the heteroclinic solutions 

cr (r) = 2arcsin(tanh(|5| 1/2 r)) + tt = 4 arctan(exp(|5| 1/2 r)), 
J (r) = ±2|5| 1/2 sech(|5| 1/2 r). 

All of these heteroclinic orbits are also heteroclinic manifolds for the corresponding unper- 
turbed stroboscopic Poincare map with period 127r/x. If 5 = 0, then the Poincare-Melnikov 
approach is not directly applicable; hence, we need to exclude this possibility in the calcu- 
lation of the Melnikov function. It is therefore clear from the results of Appendix |A] that 
all resonant circular orbits (eo = 0) are excluded from our analysis except for m = 2. This 
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is particularly significant in view of a result of linear perturbation analysis |TT[ that the 
main resonance for a circular orbit occurs for m = 2; this is consistent with the reciprocity 
between the emission and absorption of gravitational radiation. Furthermore, we need to 
exclude e$ ~ 0.76 for the (2 : 1) resonance, e$ ~ 0.85 for the (3 : 1) resonance, etc., since 
5 = at these zeros of K+(e). Let us remark here that the initial orbit is such that po > 
by convention; therefore, there is a marked difference between the absorption of positive and 
negative helicity gravitational waves. In particular, K™(e) does not vanish for < e < 1. 
Let us also note here that is an odd function of r while Jo is an even function; moreover, 
Jo — > as r — > ±oo. As mentioned previously, we will compute the Melnikov integral along 
these invariant manifolds. For regular perturbations, the existence of simple zeros for this 
function implies that, for sufficiently small //, there are chaotic invariant sets for the per- 
turbed Poincare map, and hence there are chaotic invariant sets for the perturbed system. 
In the system ( |3.13| )- fl3.14| ), the perturbation is not regular due to the rapid oscillations of 
the perturbation when ii is small. Nevertheless, we will compute the Melnikov integral for 
this system. 



4 Poincare- Melnikov Function 



Recall that for a time-periodic planar system of the form 

x = f (x, y) + fifi(x, y,t), y = g (x, y) + ngi(x, y, t), 

where the unperturbed system is Hamiltonian, if t t— > (x(t),y(t)) is a solution of the unper- 
turbed system that starts on a heteroclinic or homoclinic orbit, then the Melnikov function 
defined on this orbit is given by 

/oo 
[fo(x(t), y(t)) 9l (x(t), y(t),t + - g (x(t),y(t))h(x(t),y(t), t + £)] dt. (4.1) 
-oo 



The computation of M for the system (|3.13|) — (|3.14 ) is quite complex. However, the end 
result is reasonably simple. Let us begin with the term T = A/ Jo; as in the computation of 
formula ( |3.12| ), we have that 



_7_ 
L 



OO r 

-am K u + sin [yr] + 2g{\ )) + K v _ sin (yij - 2g(l + 



m' 



(4.2) 



Using formulas ( |3.12| ) and ( [4.2|) , let us compute 



A = L rxdg=i- l'{l-(l))dg 



such that (A) (rj) = 0; in fact, the constant of integration must vanish and the result is 



T = -am 
2 



u=l 



cos (is V + 2g(l - g 
1 - ^ 

m 



u=l 



cos [yr] — 2g(l + 
1 + ^ 

m 



(4.3) 
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Let us recall that in the definition of the Poincare-Melnikov function Q4.1|) , 

Tfi 3 <9r 

fo = Jo, fi = r - —P, g = S sinao, g x = J (-« - —), 
6 6 oa 

where 

P = — Jo + aVL{C L + -C G ), k = aVL(C lL + -C eG ). 
mm m 

Thus, we have a preliminary expression for the integrand of M (cf. Appendix |B|) 

a Q 2 dT 
fo9i ~ 9ofi = ~^Jo(CiL H Cta) ~ Jot- + 5 sin a fi- (4.4) 

The Melnikov integral ( |4.1|) is calculated in Appendix || and the result is 

M(0 = lam 2 H m I°S° 
o 

2 ~^ K v — m v + m 

ij oo 

6 ^! 

1 oo 

+ -am 2 F m J2iKe P»K + K% P~SZ). (4.5) 

Here the only quantities that depend on £ are S° = sin(£/3/x) and 

5* = sin[±(f/3^)(l T v/m) + i^r/m], 

while H m , F m , K± and K± e = dK±/de all depend upon the eccentricity of the orbit at 
resonance e = (1 — eg) a, where e = G /L . In fact, F m = (e^ — 2eo/m)/e and i? m = 
e (2e + F m ). Moreover, 1° and are also dependent upon e via 5 and are given by 



/oo / T \ 

sin ( — ) sin<T (r) dr, 
-oo v 3« y 

dr. (4.6) 



I , v 



V , s T / V 

-oo r ±— 1 T - 

m 3/i v m 



It is interesting to note that /iJ° and /i-P^ only depend upon //' = ^|5| 1/2 . A detailed 
discussion of these integrals is contained in Appendix [B|; in fact, a method is described there 
that can be used to calculate 1° and P^ whenever 2v/mis an integer. For instance, it follows 
from the results of Appendix || that P^ can be explicitly determined for m = 1,2. 

It is important to point out that the Poincare-Melnikov function is periodic with period 
6nfim and has a countable infinity of zeros at = 3/i(l + mN)ir for any integer N. 
Indeed, 5° and all vanish for ^/(3/i) = (1 + mN)ir, where N = 0, ±1, ±2, . . .; therefore, 
M(6v) = 0. 
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The zeros of the Melnikov function M(£) are all generally expected to be simple. This 
may be seen from the nature of two consecutive zeros of M(£), e.g. £o an d £1, due to the 
periodicity of the Melnikov integral. However, it is more convenient for our purposes to 
compute M'(£jv) for general N. We note that in the expression ( |4.5|) for M(£) only <S° and 

depend upon £, and 



dS° 



N) 



1 ( 

3/i 



-1) 



Nm 



(60 



T J_(_l)^(m-»/)( lT ^ 



3/i 



It follows that 



-ir m 3//M'(e 



A' J 



b d i/=i 

Y oo 

-^m^(-i)%(^ + + r?;) 

D I/=l 

+ lam 2 F m f;(-ir[(l-^ e P+ 
o ~[ m 



;i + -)^ e p-],(4.r) 

m 



assuming that term-by-term differentiation of the infinite series in ( f4.5|) is meaningful. 



Inspection of equation ( |4.7|) indicates that in general M'(^) ^ 0. To see this in some 
detail, let us focus attention on the value of M'(^ N ) for an eccentricity e such that < 
e <l. We remark that fi 2 M'(^ N ) can be regarded as a function of independent variables 
e and // = yu|<5| 1//2 ; in fact, /i' occurs only in 1° and P^. The forms of H m (e) and F m (e) 
then lead us to distinguish two cases: m / 2 and m = 2. Suppose first that m ^ 2; then, 
the discussion following the introduction of 5 ^ in equation fl3.15|) implies that e ^ 0. It 
follows from the results of Appendix [A] that the leading term in the expansion of M'(£jv) i n 
powers of eo -C 1 is 



(_l)^(m+l)( m _ 2 ) a 

6m/ie 



[(to - l)P 1 f 



TO 



3 J' 



where P+ is given by equation (|4.6|) . The general form (|B.5|) given in Appendix [B| for the 
integrals in equation ( |4.6| ) can be used to show that this leading term is nonzero for to = 1; 
a similar result is expected for to > 2. Let us next suppose that m = 2. If eo = 0, then the 
Melnikov function M(£) vanishes in this case; further analysis of this limiting case is beyond 
the scope of this paper. For eo 7^ 0, the results of Appendix ^ imply that the leading term 
in the expansion of M'(£jv) m powers of eo <C 1 is now 



l) N ae 
6/i 



25P 3 + ), 



where P^ and P 3 + are given by equation ( |4.6| ) for m = 2. The general form ( |B.5| ) given 
in Appendix [B| for the integrals in equation (}4.6|) again indicates that this leading term is 
nonzero. Hence, M(£) has in general simple zeros for < eo <C 1. 
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5 Numerical Experiments 



We now consider the same Hamiltonian as in the calculation of the Poincare-Melnikov func- 
tion, 

H = l(Pr + $) ~ 7 " + COS 20, (5.1) 

and will demonstrate numerically that this system could be chaotic near a resonance. We 
first fix an energy surface H(p r ,pg,r,9) = h and consider all orbits confined to this energy 
surface. For instance, choosing an orbit with (p r ,pe, r, 9) = (e, 1, 1, 0) at t — fixes 

The equations of motion are 

dr 

~dl 
d6 

~dl 

dp r 

~dt 

^ = eaQ 2 r 2 sm29, (5.2) 
dt 

which can be integrated with different initial conditions all with H(p r ,pg,r,9) = h. This 
energy equation can be used to calculate pe algebraically in terms of the other variables; we 
limit our attention here to the solutions of this quadratic equation with pe > 0. Moreover, 
we choose 9 as our independent variable and note that the resulting system for (p r ,r) is 
periodic in 9 with period it. We then consider the Poincare map and plot in the (p r , r)- 
plane the result of our numerical integration of the system at 9 = 0, ir, 2ir, .... In following 
this procedure, we encounter the difficulty that 9 is not necessarily monotonic. Despite this 
complication, it is possible to obtain useful information from the integration of system ( f>.2\ ) 
as in our previous work (see figure 1 of To avoid this problem altogether, let us restrict 
attention to those orbits that stay with a simple branch of pg > 0. That is, we note that 
H(p r ,p 9 ,r,9) = h has the solution 

p e = lnr 2 [l±(l-U) 1 /% 



Pr, 




Pe ^ 




r 2 2' 




P±_h_ 


eaVL 




eafl 2 r 2 sin 


29, 



where 



U=^4p 2 r-y-h)+4eacos29. 



Thus the system (|5.2j ) reduces to 



~dt =V - 



dr 
It 

% = -IpJ + ^ + ^ + Itf r [l - 4m cos 29 ± (1 - Uf\ (5.3) 
dt r r z r 2 



15 



where the upper (lower) sign corresponds to the positive (negative) branch of 9. We have 
subjected the system (|5.3|), restricted to the positive branch of 9, to numerical analysis and 
have found numerical evidence for Hamiltonian chaos as presented in figure 2. 



It is interesting to consider scale transformations in the dynamical systems ( |5.2|) and ( |5.3|) 
along the lines described in detail in our previous work || [|, |TI]] . The result is that we can 
use dimensionless quantities and set ko = 1 once we use the initial semimajor axis of the 
binary system as our unit of length and the initial period of the binary divided by 2tc as 
our unit of time. We shall use this convention in the numerical experiments reported in this 
paper. 



Partial phase portraits of the Poincare map for system (|5.3j ) are depicted in figures 2 
and 3, where we have chosen the parameter values (k = 1) 

n = 1, a = 2, (5.4) 

for the plane right circularly polarized gravitational wave that propagates perpendicularly 
to the orbital plane. Also, for the top panel of figure 2 the external perturbation is absent 
so that 

e = 0, h = -0.999982, (5.5) 
while for the bottom panel of figure 2 and for figure 3 

e = 0.0131, h = -0.986882. (5.6) 

In these graphs, the horizontal axis corresponds to the variable p r while the vertical axis 
corresponds to r. Inspection of the system (|5.3| ) reveals that the phase portrait of the 
Poincare map is in general symmetric about the line p r = 0. 

The top panel of figure 2 depicts several orbits corresponding to invariant tori for the 
unperturbed system (e = 0), whereas the bottom panel depicts several orbits for the per- 
turbed system (e ^ 0). Let us note that the stochastic region in figure 2 is obtained as a 
single orbit. This same orbit is also shown in the top panel of figure 3. It can be obtained 
via integration of system ( ^3|) , or (|5.2| ), with 



(r,0) = (0.61927711963654,0) 

and 

(p r ,p e ) = (0.51405620574951,0.83454334164660) 



as the initial conditions and with the parameters given in equations (|5.4j ) and (|5.6| ). This 



chaotic orbit is near the unperturbed (1 : 1) resonant torus. Let us recall here that the 
energy of the osculating ellipse is given by 

so that the semimajor axis is a = —ko/(2E) and the corresponding Delaunay element is 
L = {ha) 1 ' 2 . The ( m : n) resonance occurs at mk^/L^ = nQ; however, in the case under 
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Figure 2: The phase portrait of the Poincare map for system ( |5.3|) with the parameters given 
in displays ( |5.4| )- (|5lj| ). The top panel is a phase portrait for the unperturbed system ( |5.5|) 
while the bottom panel depicts a phase portrait for a perturbed system (|5\6| ) on a nearby 
energy surface. 
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Figure 3: The top panel shows a single orbit of the Poincare map for system ( |5.3| ) with the 
parameters given in displays ( |5.4|) and fl5.6|) . The bottom panel is a graph of the Delaunay 
variable L versus time for the same orbit. It follows from this plot that the behavior of the 
semimajor axis (a = L 2 ) of the osculating ellipse is chaotic in this case. 



18 



consideration here the primary resonances are (m : 1) resonances and our chaotic orbit is 
near the dominant (1:1) resonance. Thus L for the chaotic orbit is likely to fluctuate about 
the (1 : 1) resonant value Lq = 1, though in the case under consideration we expect that 
secondary (m : n) resonances with n > 1 are also involved since e is large enough. The 
bottom panel of figure 3 is a graph of the Delaunay variable L versus time for the same orbit 
depicted in the top panel. The stochastic region filled out by the orbit depicted in figure 3, 
an orbit of a 2-degree-of-freedom time-independent Hamiltonian, seems to be bounded by 
KAM tori. The fine structure of the graph of L versus t over short time scales, i.e. a few 
thousand units of time, shows both chaotic and apparently regular oscillations. These are 
accounted for by reference to the orbit of the corresponding Poincare map shown in the 
top panel of figure 3. The corresponding trajectory visits regions near hyperbolic saddle 
points where transversal crossings of stable and unstable manifolds induce a homoclinic 
tangle with embedded Smale horseshoes, etc. Thus a sojourn near such a region produces a 
highly chaotic apparent change in the semimajor axis of the corresponding osculating ellipse. 
However, when the orbit is "distant" from these hyperbolic saddle points, the semimajor axis 
of the osculating ellipse oscillates with apparent regularity though the regimes of apparent 
regularity are dependent on the scale at which the fluctuations in L are observed. This 
results in a form of Hamiltonian intermittency. Of course, the sojourn time for apparently 
chaotic (large-amplitude chaotic signal) versus apparently regular (small-amplitude chaotic 
signal) oscillations is itself highly chaotic in our numerical experiments. 

The width of the stochastic layer (SL) in figure 3 may be estimated on the basis of the 
following considerations: If e = 0, one finds that the variation of E — h + Qp$/2 vanishes 
along the orbit (cf. eq. Q5.2| )) and hence L is constant (5L/L = 0). Physically, this comes 
about since an ellipse is simply viewed from a rotating frame. Once the perturbation is 
turned on, the KAM theory implies that for sufficiently small e > the motion is bounded 
for all time. Inspection of the system ( |5.3| ) then reveals that the small parameter in the 
problem is effectively 4\a\e as < | cos2#| < 1; therefore, we expect that 5L/L « 4|a|e in 
agreement with the data of figure 3. 

The phase portrait of the Poincare map for system ( |5.3| ) exhibits extreme sensitivity to 
changes in the parameters of the system. For example, stochastic regions comparable in 
relative size to the one depicted in figure 3 are not easily detectable for a random choice of 
parameter values. We have verified all of the main numerical results reported in this paper 
via two different stiff integration routines using double precision arithmetic. We emphasize 
that a detailed numerical study of the Hill system is beyond the scope of this paper, since 
the relevant parameter space is rather large. Our numerical results nevertheless support the 
view that there is chaos in the Hill system. On the other hand, it would be interesting from 
an astrophysical standpoint to see the chaotic behavior of the binary in the inertial frame of 
reference. This issue is taken up in the next section. 
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6 Inertial Observers and the Hill System 



The equation of motion (1) can be derived from the Hamiltonian 

U = -P 2 - ^ + -KaX^ 
2 r 2 3 

in the inertial frame of reference. This Hamiltonian can be written in polar coordinates 
(r, 9) as 

H = \{ P r + pF) - 7 + l r2 tt° + ^+ COS ( 20 " «*) + £~ Sin ( 2@ " (6- 1 ) 

where we have used equation (2) for the tidal matrix. Let us recall here that Q can be 
positive or negative. 

The connection between the inertial motion and the Hill motion based on Hamilto- 
nian ( |2.6| ) can be simply obtained by referring the inertial motion based on the Hamiltonian 
(|6.1|) to the rotating frame. In this process, = 6 + Qt/2 and the canonical momenta remain 
invariant, i.e. P r = p r and Pe — Pe, while the Hamiltonian in the inertial frame is related to 
that in the rotating frame by the well-known relation 

H = H+hn Pe (6.2) 

that expresses angular momentum-rotation coupling. It is important to observe that with 
the initial conditions (p r ,pg,r,9) at t = being common to both the inertial and rotating 
frames, we can integrate the equations of motion in these frames starting from the same 
initial osculating ellipse. 

Let us now choose £o — £- = and £ + = eaQ 2 as before. It then follows that the 
perturbation in the inertial Hamiltonian is of the form |eafi 2 r 2 cos(26 =F &t), where the 
upper (lower) sign stands for plane right (left) circularly polarized gravitational radiation of 
frequency Q > that is normally incident on the Keplerian binary system. To describe the 
inertial motion, we introduce the osculating ellipse with semimajor axis a, eccentricity e, 
etc., in the inertial frame and the associated Delaunay elements (L,G, £,g). The equations 
of motion in terms of Delaunay's variables are 

dL _ c &H ext 
dt dl ' 

^5 = c dH ext 
dt dg 
dl = k 2 &H ext 
dt L 3 dL ' 

f= e ^fet, (6.3) 
dt dG K ' 



and 



— = a 
dt ' 
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where r = f2t is a new variable and 

H ext (L,G,£,~g,f) = ^Vcos(2e T f). 

We are interested in the general behavior of system (|6.3| ) taking into consideration the fact 
that the motion is always bounded for e < ej <;am as a consequence of the KAM Theorem. 

There are two frequencies in the dynamical system under consideration here, namely 
the Keplerian frequency uo = k^/L 3 and the frequency of the external radiation Q; hence, 
the possibility of resonance cannot be ignored. Thus the averaging principle is in general 
inapplicable here due to the fact that resonance could occur. Indeed, resonance could come 
about when uj and Q are commensurate, i.e. relatively prime integers m and n exist such 
that muj = nQ; at this (m : n) resonance, L is fixed at its resonant value Lq. It will be shown 
in this section that for system (|6.3| ) the primary resonances are in fact (m : 1) resonances, 
i.e. n — 1. The resonant behavior of the system around the (m : 1) resonance manifold can 
then be described by the partially averaged equations for system (|6.3p . 

It is important to point out here that the (m : n) resonance is invariant, i.e. the resonance 
is the same whether it is encountered in the rotating frame of reference or in the inertial 
frame of reference. To make this point explicit, let us note that the energy of the osculating 
ellipse 

2 r r 2 r 

is an invariant quantity in the sense described above {£ = E), and hence so is the semimajor 
axis of the resonant orbit and its corresponding Delaunay element L = L = (— /cq/2^) 1 / 2 that 
becomes fixed at resonance. The same holds for the angular momentum (G = G) and the 
eccentricity (e = e) of the osculating ellipse; nevertheless, the tildes on the Delaunay action 
variables will be generally kept in the rest of this section as a reminder that the Hill system 
is being viewed from the inertial reference frame. Chaos is invariant as well and expected 
to occur near a resonance; indeed, this point can be illustrated with a single chaotic orbit 
from the case studied in figure 2. This chaotic orbit that is near a (1 : 1) resonance is given 
in figure 3. Further illustration of small-amplitude chaos is contained in figures 4 and 5 
described below. Figure 6 illustrates large-amplitude chaos that comes about when primary 
resonances overlap. 

The average behavior of the system around the (m : 1) resonance manifold can be 
obtained from the second-order partially averaged dynamics. To this end, let 

A) 

such that the second-order averaged equations in the inertial frame are given by 0, |l(J 

V = —ez(—mT c sinm<l> + mT s cosm<3>) 

^/ dT c . ^ dT s 
— eD[—m — =- sinm<P + m — — cosmffl , 
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~ dT BT 

G = — el -— cos m<P + -— sm m<P ) , 
dg dg 

■ i3k 2 V ^klV 2 dT c dT s 

$ = -ea^- h e( — Sr- 1 cosm$ H =- smm$), 

Lq L 5 8L dL 

x. dT c dT s 

q = e( — ^cosm$H ^-sinm$). (6.4) 

dG dG 

Here T c = J^(L, G) cos 2g and T s = T^±(L, G) sin 2g, where 

jrm = ^Ltfa 2 K™{e) (6.5) 
2m 

for n — 1. If n ^ 1, T c = T s — and primary resonances do not occur. It follows that for 
system ( |6.3| ) the primary resonances are given by mu = Q. The system ( |6.4|) is evaluated at 
L = L . 

The second-order averaged dynamics can be simplified if we introduce a new variable 
A = m$ ± 2g. It is then straightforward to recast ( |6.4[ ) into the form 

V = me l (jr™ + e^P^ £ ) sin A, 
G = ±2e^sinA, 



A = - e s + e—V' z + e{mF™ L ± 2JF^) cos A, (6.6) 



where JF™^ = dJ 7 ™ /dL, etc. Restricting attention to the first-order averaged equations in 
( |6.6|) , we find that this system can be integrated. Thus 

G-G Q = ±-(L-L ), (6.7) 
m 

where Go is the orbital angular momentum when the system is exactly at resonance. More- 
over, 

CvTT) — 

V 2 = — Lg#£(eo)(cos A - cos A ), (6.8) 

where Ao and eo belong to the relative orbit at resonance. It follows from equation ( |6.8|) that 
T> has oscillatory character and the maximum value of T> 2 occurs at either A = (0, 2n, 4n, . . .) 
or (7r, 3n, 57r, . . .) depending upon whether aK±(e ) is positive or negative, respectively. The 
amplitude of the total variation in L is thus given by SL = 2e 1 / 2 r , max, where 

^rnax = (e~o) sin 2 ^% ± g ) (6.9) 

for aK± > 0, and 

^max = -— L 2 K^(e Q ) cos 2 (-$ ± g ) (6.10) 

for aK™ < 0. Other properties of the motion can be studied as well since A varies just like 
a standard nonlinear pendulum. In addition to these results from the first-order averaged 
equations, small regular corrections exist that are beyond the scope of this paper and are 
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due to terms of second order in e 1//2 in the averaged equations. Furthermore, certain chaotic 
effects are also expected near the primary (m : 1) resonance. 

These theoretical conclusions may be illustrated by a simple example: Imagine a resonant 
orbit with (p r ,Pe, r , 9) = (e, 1, 1, 0) at t = such that e = 1/2, g = —n/2, u = tt/3, v = n/2 
and hence $o = tt/3 — 3^ 2 /4 for this orbit. We choose k = 1, e = 1CT 3 , a = 2 and 
fi = m/Lp = m(3/4) 3//2 , so that the orbit is initially in (m : 1) resonance. The results 
of the numerical integration of this system for m = 1,2 are in good agreement with the 
predictions of the averaged equations (|6.7|) and ( |6.8| ), since the amplitude of the chaotic 
motion that is present near the resonance is smaller than, or at most comparable with, the 
amplitude of resonant motion. This is illustrated for the (2 : 1) resonance in figures 4 and 5 
for the incident right and left circularly polarized (RCP and LCP) waves, respectively. It is 
clear from these results that there is a high-frequency component with dominant frequency 
« Q that is at least partly chaotic and is superposed on an average regular motion. The 
net amplitude of this "fast" component can be crudely estimated to be 5L/L w 0.005 and 
5G ~ 2 5L regardless of polarization; this is in rough agreement with figures 4 and 5. For 
the regular motion A /2 = $ ± g ; hence, for the RCP case A /2 = — (n/6 + 3 1/,2 /4) and 
if! (0.5) ~ 0.923 so that using equation ( |6.9| ) we find 5L « 0.094 in good agreement with the 
numerical results in figure 4 when the small-amplitude chaotic signal is ignored. Similarly, 
for the LCP case A /2 = 5vr/6-3 1 / 2 /4, (0.5) « -0.009, and using equation ( |6U0|) we find 
5L « 0.006 for the amplitude of regular oscillations as in figure 5. It is interesting to note 
that the amplitude of the "chaotic" signal in L is about the same for RCP and LCP waves. 
However, the response of the orbit with pg > to LCP radiation is smaller by an order 
of magnitude than its response to RCP radiation; hence, the "chaotic" signal in figure 5 is 
comparable in amplitude to the regular signal that is consistent with the first-order averaged 
equations. Thus the chaos that appears in figures 4 and 5 has small amplitude; in fact, to 
encounter large- amplitude chaos as in figures 2 and 3, it seems from our numerical work 
that the strength of the perturbation must be larger by at least an order of magnitude. The 
variation of the angular momentum G is equal to the variation of L for m = 2 and in (out 
of) phase with it for RCP (LCP) waves in accordance with equation ( |6.7|) and in agreement 
with the numerical data of figures 4 and 5. Furthermore, it is clear from equation ( |6.6j ) that 
the angle A satisfies the standard equation for a nonlinear pendulum to first order in e 1 / 2 
and the frequency associated with the small-amplitude oscillations of this pendulum is Co 
given by 

u 2 = ^emQ 2 \aK^(e )\. (6.11) 

For the RCP case depicted in figure 4, we find that 2tt/uj « 65. The pendulum involved here 
is nonlinear with amplitude Ao(RCP); therefore, the relevant period of oscillation is longer 
by about 29% resulting in a predicted period of m 84 in agreement with the period of regular 
oscillations depicted in figure 4. Similarly, for the LCP case depicted in figure 5, we find 
from equation ( p. 11] ) that the period of small-amplitude oscillations is ~ 670. The period 
of regular oscillations can again be calculated for the nonlinear pendulum with amplitude 
A (LCP), and it turns out to be longer by about 10% in this case due to the nonlinearity. 
Thus the predicted period is ~ 739, which agrees with the numerical results based on the 
exact system given in figure 5. 
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Figure 4: The behavior of relative orbit in a binary system with initial conditions 
(p r ,pg,r,6) = (0.5,1,1,0) in (2 : 1) resonance with a normally incident right circularly 
polarized (RCP) gravitational wave. Here k = 1, ea = 0.002 and = 2(3/4) 3 / 2 « 1.299. 
The top panel depicts L versus time, the middle panel depicts the angular momentum G 
versus time and the bottom panel depicts the eccentricity e versus time. 
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Figure 5: Same as in figure 4, except that the incident radiation is left circularly polarized 
(LCP). 
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Large- amplitude chaos near primary resonances can still occur in our system if the 
Chirikov resonance-overlap condition holds. For instance, numerical experiments with the 
initial orbit under consideration here indicate the presence of large-scale chaos if the incident 
RCP wave has frequency Q = tc/2. To explain the appearance of large-scale chaos in terms 
of the Chirikov criterion, let us note that a primary (m : 1) resonance would correspond to 
a Delaunay variable L m in this case given by L m = (2m/7r) 1 / 3 . Thus the initial Delaunay 
variable L of the orbit in our example, L = (4/3) ^ ~ 1.15, falls between those of (2 : 1) and 
(3 : 1) resonances since L 2 ~ 1.08 and L 3 « 1.24. The Chirikov criterion for large-amplitude 



chaos [14], [Tj| can be written in this 

^1/2 _ L m +i 



ch [(2m/3)P m \aK^(e Q W^ + [(2(m + l)/S)Ll +1 \aKr\e )\}^ ' 

using equations Q and flOpp . From K\{Q.h) « 0.92 and K% (0.5) « 0.76, we find that for 
m = 2, 

e ch « 1.65 x 10~ 3 . 

Our numerical computations suggest, as is expected, that large-scale chaos occurs for a 
somewhat smaller choice of the perturbation parameter, certainly for e = 10~ 3 . At much 
smaller e, for example e = 10 -5 , the action L still exhibits large-amplitude chaotic effects, 
but seems to be bounded by KAM tori. However, for e > it appears that the action L 
is no longer confined by KAM tori. This is illustrated in figure 6, where the nature of the 
binary orbit is studied for Q = n/2 and e = 10~ 3 . At this Chirikov threshold, all KAM tori 
are expected to disappear and gravitational ionization might take place. It should be noted 
that in our previous work ||, the ionization process was associated with Arnold diffusion. 
Here, however, gravitational ionization is due to the fact that e > e^anr In figure 6, the 
perturbed orbit is displayed by plotting Y = rsinG versus X = rcosO. In the absence of 
the perturbation, r varies from a(l — e) = 2/3 to a(l + e) = 2 along the initial elliptical orbit; 
however, r can reach relatively large values after the RCP perturbation is turned on. In the 
experiment depicted in figure 6, for instance, the orbit stays generally close to the initial 
osculating ellipse for about 100000 time units, but large-scale deviations set in eventually 
such that r can have values as large as ~ 135 over an interval containing approximately 
450000 time units. 

As in our previous work on gravitational ionization [T], [7], |8| , we find that the binary system 
is not completely destroyed in the process of gravitational ionization; that is, one member 
of the binary does not go away once and for all. Rather, the system alternates between 
dissociation and recombination as in figure 6. This is essentially due to the oscillatory 
character of the energy exchange between the wave and the relative orbit. In a realistic 
astrophysical environment, however, new phenomena might lead to the complete ionization 
of the binary system once one member ventures sufficiently far from the other one. 



Finally, let us assume that the system ( |6.3| ) is off resonance, i.e. it is far from the primary 
(m : 1) resonance manifold. Then we expect that the motion would involve small amplitude 
oscillations of frequency ~ Q and amplitude ~ e near the initial osculating ellipse. We have 
explicitly verified this point for the example given above but with Q = 7r/10. 
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Figure 6: Plot of the perturbed orbit exhibiting large-amplitude chaos. The initial conditions 
are (p r ,Pe, f, 9) = (0.5, 1, 1, 0) with k = 1, and the parameters for a normally incident RCP 
wave are given by e = 10~ 3 , a = 2 and Q = tt/2. The top panel depicts a chaotic rosette 
corresponding to the relative orbit of the binary system over approximately 450000 time 
units, the middle panel depicts the energy £ = E given by equation (|5.7p versus time and 
the bottom panel depicts the orbital angular momentum pe = G versus time. 
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7 Discussion 



The KAM analysis of Hamiltonian chaos is explicitly illustrated in this paper for the in- 
teraction of normally incident circularly polarized gravitational radiation with a Keplerian 
binary system. In the rotating Hill coordinates, the gravitational wave stands completely 
still and hence the Hill system becomes autonomous. The existence of chaos in this system 
is demonstrated numerically. On the analytic side, there is no proof at present that the stan- 
dard approach to Hamiltonian chaos — i.e. the Poincare-Melnikov theory — is applicable 
to the Hill system. Nevertheless, the Melnikov function has been useful in this regard in 



certain analogous circumstances [12], |13[. Therefore, we compute the Melnikov function for 
the Hill system and show that it has the requisite properties for chaos. The astrophysical 
implications of these results are explored; in particular, we employ the partially averaged 
equations to describe the nature of the relative orbit when the binary is at resonance with 
the incident wave. The possibility of gravitational ionization in the Hill system is discussed 
in connection with large-scale chaotic orbital motions brought about by incident radiation 
with an amplitude e > e c ^, where corresponds to Chirikov's resonance-overlap criterion. 



A Some Properties of K±(e 

The quantities 



K™(e):=^m(A m ±B m ), m = 1,2,3,..., 



can be expressed in terms of Bessel functions via equations ( [2.1 0|) and ( p.ll| ). Using standard 
expressions for the Bessel functions, we find that for e C 1, 

Thus, K\{e) = -3e + 0(e 3 ), K%(e) = 2 - 5e 2 + 0(e 4 ), K\{e) = 3e + 0(e 3 ), etc., for e -> 0, 
so that K™(0) = 25 m2 . Moreover, for e < 1, 

k - w = - ( ^SS em+ ' + ° (em+1)i (A - 2) 

hence, K™(0) = 0. It follows that K™(e) < on the interval < e < 1, since we have 
already shown in [jl| that K™(e) is monotonically decreasing on this interval (see p. 107 
of i). 

We are also interested in the function K™ e = dK± / de. It is simple to see from ( |A.1|) that 
K\ e (0) = -3, ^ e (0) = 3, and for all other values of m ^ 1,3, K™ e (0) = 0. Similarly, it 
follows from ( |Q|) that K™ e (0) = 0. 

For e — > 1, one can use the definition of K± to find 

K%(1) = --J m (m). (A.3) 
m 



28 



Here Ji(l) « 0.44, J 2 (2) « 0.35, etc., and for m > 1 

~ -T§73y m " 4/3 - ^ 

It is also interesting to determine the behavior of K± e as e — > 1. We find that K± e (l) 
diverges as 

KZ(e) „ ±(l-e 2 )- 1 / 2 ^J' m (m), (A.5) 

e— >1 

which agrees with the results presented in figure 1. 



B Calculation of M(£ 

Let us consider the solution r i— > ( Jo(r), co(r)) for the unperturbed heteroclinic orbit of 
the system ( |3.15|) . Starting with equation Q4.4|) , we note that the function r i— > sin cr (r) 
is odd, so only the odd terms in the expression for fx contribute to the Melnikov integral. 
Likewise, the function r i— > Jq(t) is even, so only the even terms of Cn + (2/m)Cic an d 
dT /da contribute. To compute (3 in the formula for fx, let us consider C in the form 



5 L 4 / G 2 



C(L,GJ,g) = -^{l-—) cos 2g 



T 4 CO 1 



and determine its partial derivatives with respect to the action variables 



C L (L, G, £, g) = ^{2L 2 - G 2 ) cos 2g 



3 oo 



LI/ 

+ 4-ra £ ~ cos(2c/ + z/£) + iC cos(2c/ - I/€ 



T 4 OO 1 

+ T2 E - cos ( 2 2 + v£) + cos(2^ - 

K u=l V 



and 



C G (L,G,£,g) = -5—cos2g 



I?G 



r 4 co i 

+ TT E - cos ( 2 # + ^) + #-g cos ( 2 9 - vi 

^0 i/=l v 

where K+ L = dK+/dL, etc. Also, let us note that here £ = (ao + n — 2g)/m and g = r/(6/i). 
Of course, we must replace r by r + £ in the Melnikov integrand wherever there is explicit 
dependence upon time in fx and g%. For notational convenience, let us define 

a = ^(r)±ifi=F-). 



m 3/i v nv 



S± = ±f it- + 
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We then have an expression for fi = T — (m/6)/3 using equation 

1 ~ C0s(C+ + 5 ± ) 1 ^ COS(C-+0 



5 „, IGo lGgx , T £ 
-am (1 --y) cos ( — + — 



2 OO 1 

-am 2 y - [Kl cos(C+ + s+) + cos(C- + s. 

3 f V 



1 00 1 

- -am 2 L J2 ~ (Kl cos(C+ + s+) + fC L cos(C_ + s. 
1 00 1 

- -amL y ~(K+g cos(C+ + s+) + #-g cos(C- + s_ 
^ i/=i ^ 

Moreover, the odd part of f\ is given by 

x odd 1 t^u sin C+ sin s + 1 ^ sin C_ sin s_ 
f, u = — am > iTf — + -am > K 

u=i rn v=l m 

+ -am (1 -J j sin (— ) sin (— ) 

3 v m L 2Lq' 3/i 3/x 

2 00 1 

+ -am 2 ^ sm C+ srn s + + si n C- srn s -) 

3 „=i v 

1 00 1 

H — am 2 L —(KK L sin£ + sin s + + K V _ L sin sin s_) 
6 ^ 

1 00 1 

H — amL ^ —(K+ G sin £ + sin s + + i£^ G sin £_ sins_). 

3 „=1 ^ 

On the other hand, only the even terms of Cn + {2/ m )CeG contribute to M; to compute 
these terms, we differentiate the expressions for C L and C G with respect to £, 

AT3 oo 

C tL = -72" E ( - ^+ Sill ( 2 ^ + + Sm ( 2 # - ^)) 
r 4 oo 

+ T2 E ( " K Il sm(2^ + ^) + K V _ L sin(2^ - W 

T 4 oo 

= T2 E ( - ^G sm(2^ + ^) + JC G sin(2 5 - z/£)) . 

fc „=1 

Thus, we find that at resonance 

C tt + — CtG = -77- E sin (C+ + s+) + K v _ sin(C_ + s. 



mL cx 



n 



£ sin(C + + «+) + sin(C_ + «_ 



i/=i 





2L 00 



n ,=1 



^ fc G sin(C + + 8+) + K»_ G sin(C- + s_ 
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and the even terms are 



(Cil H C-ig) 

• m ' 



4m 00 

-—— ^2 {K+ cos C+ srn s + + K- cos C- si n s -) 
^ i/=i 



y)(K+L cos C+ si n s + + K-l cos C- srn s -) 



00 



2J (-K+g cos C+ si n s + + K-g cos C- srn s - 



Finally, we must compute — Jq dT/da, where 



dr 

da~ 



-a 



£ ^ sin(c+ + g+) - f; vki sin(c - 1 s 



1 + = 



using equation 
integral, hence 



Only the even terms in dT/da would contribute to the Melnikov 



/<9I\even 1 



„ cos C+ sin s_ 



i/=i 



cos£_ sins. 

m 



We now have 



1 f 00 ^^1 I r°° / 2 \ even 

M(0 = --amK™ sin a (r) f° dd dr + -atl J^(r)(C eL + -C eG ) dr 

O J— 00 O J— 00 v 771 7 



Let us define 



and note that 
It follows that 



2 /driven 



— amK' 1 
6 



sino"o(r) sin f— 1 dr, 



r± 



M(0 



1 r 00 
--amK™ / sin o"o(r) sin o?r, 
6 J —00 

5 :=sm— , S„ := sms±, 
3/z 

/OO 
Jq(t) cos C±^t, 
-00 

r+ — _ r~ 7+ — r~ <?+ — _ 



•-am i— 7,t«S,t H — am = —I~S~ 

— "=l m ~" i/=l ~ m 
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3 v m L 2 Lq > 

2 00 1 

6 v=\ V 

1 00 1 

1 00 1 
+ g«^o E T^nSt + K^I-S-) 



V 



cy OO 

1 OO 

1 OO 

- -aL Y,(K G J+S+ + K^ G r v S-) 
1 00 vK v 1 00 i/if" 



(B.l) 



Furthermore, let us define 



<9e d 



«9G 
<9e 



K±L ~ dLde K± - dL ±e ' 



F m (e) := L 



de 2 de 
dL m dG 



(Lq,Go) 



We recall that e 2 = 1 — G 2 /L 2 ; hence, 



de G 2 de 



G 



6 dL L 3 ' 6 dG L 2 ' 



so that 



and thus 



de_ J_ 
d~L~~eL 



F m (e) 



(1-e 2 ), ^ = -^(1-^)^, 



<9G eL 



1-e 2 --(1-e 2 ) 1 / 2 
m 



Moreover, it is useful to define 



H m (e) :=l+e 2 --(l-e 2 ) 1 / 2 , 



m 



so that 



# m (e) - eF m (e) = 2e 2 . 
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Finally, let us define 

P ± — -I* - — J± 



V 

v m 



After collecting terms in equation ( |B.1| ) and using the definitions given above, we have an 
expression for the Melnikov function: 

M(0 = lam 2 H m (e)I°S° 
o 

1 ~ uKl 1 ^ vK v 



-am y ^P+St + -am V — —- P~S. 

u=l m v=Y m 



_|_ v_ V V 



r) OO 

+-am 2 Y,(KKs+ + kip-s;) 

1 oo 

+ -am 2 F m (e) £(i^ e P+<S+ + K^P'S'). (B.2) 
t/=i 

Let us note that in this expression s + (u = m) = n so that = 0. Furthermore, we have 
from the definition of 1^ and the fact that 5sino"o(r) = —dJo/dr, 



h = / -— smC±ar= / J c 



oo 

? 



m 3/i v m 



cos (± dr, 



via integration by parts (Jo — > as r — > ±oo). It thus follows that P^ = v~ l I„ — m _1 J^ 
is given by equation ( |4.6| ). It is interesting to note that in the formula for M(£), = 0; 
moreover, the quantity 

-± / J (r)cosC±^r (B.3) 



1 =f i//m 3^/x 

is well defined even for v = m. In fact, we find from equation ( |3.15| ) that 

p+ l ,00 

lim —. — = / Jo(f) cos do rfr = 0. 

v->m i _ v j m 3772^ 7„ 00 

Using these results in equation (|B.2|) , we recover equation (|4.5|) for M(£). 

The function M(£) refers to a resonance between a Keplerian orbit and an incident 
circularly polarized wave. Therefore, M(£) depends upon the wave amplitude a, the order 
of the resonance m, the orbital eccentricity at resonance eo and the perturbation parameter 
//. The explicit form of M(£) involves certain combinations of the Bessel functions (i.e. 
K± and i£± e ) that depend only on the eccentricity and the integrals in 1° and P^ . The 
complete evaluation of P^ is beyond the scope of this work; however, P^ can be computed 
using contour integration when 2u/m is an integer. To see this, let us assume that 5 = 
— (am/6)K™(eo) > and choose the principal branch of (Tq(t) in the following. Hence the 
integral in P^ takes the form 

/oo 
sech (x) cos[j arcsin(tanha;) + wx] dx, (B.4) 
-oo 
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where j = 2v/m and w = ±(1 =F z//m)/(3/i<5 1 / 2 ). For the analytic extension of the integrand 
to the complex plane, it is essential to rewrite the integral in the form 

(sechx) J+1 [(l + i sinhx)- 7 exp(iwx) + (1 — i sinhx)- 3 exp(— iwx)} dx, 

-oo 

that is appropriate for the principal branch of a under consideration here. We note that 
the integrand is an even function. Let us now imagine a rectangular contour in the complex 
(x, y)-plane that is symmetric about the y-axis with the two long (eventually infinite) sides 
parallel to the x-axis at y = and y = 2n. The singularities of the integrand within this 
contour occur at in /2 and 3m /2. These are poles once j is a positive integer. In this case, 
a standard application of the Cauchy Theorem for this contour results in 

[1 — cosh(27ru7)] Xj W = 2m [residue (in/2) + residue (3i7r/2)], 

where the symmetry of the integrand has been taken into account. The calculation of the 
residues for arbitrary j is straightforward, but will not be attempted here; however, it is 
possible to show that Ij W has the general form 

_ exp(— ixw/2) . , . , 

T jw = 4tt W y \. \ V > , B.5 
smh(7rw) 

where Wj(w) is a polynomial in w of degree j '• — 1 such that jWj(0) vanishes for even j 
and is equal to (— l)^'" 1 )/ 2 for odd j. We find that Wi(w) = 1, W 2 (w) = —w, W 3 (w) = 
(— 1 + 2w 2 )/3, etc. Note that Ij W is well defined for w — > 0. Moreover, J° can also be 
computed using the approach outlined above. We find that 

■i \ ■ i \ a 2nWo 
sm^o^J sma"o(a:J ax 



cosh(7rwo/2) 

for the principal branch of er . For J°, w = l/(3fi5 1 ^ 2 ); hence, 

" 3 / ucosh[7r/(6^ 1 /2)]- ^ - u > 
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